System and method for modeling shale discontinuity

ABSTRACT

A method and system of locating discontinuities in a reservoir are provided. The method includes identifying a plurality of injection and production sites in the reservoir, injecting fluid into selected injection sites in accordance with an injection schedule, monitoring an output at the production sites, determining a time delay between the injection sites and the production sites based on the monitored output, and determining a location of discontinuities based on the determined time delay.

CROSS-REFERENCE TO RELATED APPLICATIONS

The present application is based on and claims priority to U.S.Provisional Patent Application No. 61/586,370, filed on Jan. 13, 2012.The entire content of which is incorporated herein by reference.

BACKGROUND

1. Field

The present invention relates generally to modeling reservoircharacteristics and more particularly to detection or prediction ofshale discontinuities based on reservoir flow data.

2. Background

A good geological model is critical to model the fluid flow in awaterflood project. Especially when we have limited water resources, thechoices of injecting locations can become an important issue inwaterflood management. For a layered reservoir with interlayer crossflow, the earlier depletion of the higher permeability layer ispredictable. In this situation, it may be useful to change the injectionschedule to enhance the oil recovery in the lower permeability layer. Toachieve this goal, identifying the shale discontinuities can be helpfulin modeling cross layer fluid flow.

In the past few decades, many different methods have been developed toestimate the parameters for layered reservoir. Seismic crossholetomography (Chapman and Pratt, 1992) can provide very high resolutionresults if enough sensors are installed. Tracer tests (Vasco et al.1999) are also widely used to investigate the flow property in thereservoir.

However, all of the above methods require additional cost and mayinterrupt the daily operation. Seismic crosshole testing usually takes along period of time and has difficulty in characterizing the flowproperty directly. Tracer tests can map the inter-well flow propertyconveniently, but usually are not repeatable.

SUMMARY

An aspect of an embodiment of the present invention a method of locatingdiscontinuities in a reservoir. The method includes identifying aplurality of injection and production sites in the reservoir; injectingfluid into selected injection sites in accordance with an injectionschedule; monitoring an output at the production sites; determining atime delay between the selected injection sites and the production sitesbased on the monitored output; and determining a location ofdiscontinuities based on the determined time delay.

An aspect of an embodiment of the present invention includes a systemfor locating discontinuities in a reservoir. The system includes aprocessor configured to: (a) determine a time delay between injectionsites and production sites in a the reservoir based on monitored outputat the production sites from injecting fluid to selected injection sitesin accordance with an injection schedule; and (b) determine a locationof discontinuities based on the determined time delay.

Other aspects of embodiments of the present invention include computerreadable media encoded with computer executable instructions forperforming any of the foregoing methods and/or for controlling any ofthe foregoing systems.

BRIEF DESCRIPTION OF THE DRAWINGS

Other features described herein will be more readily apparent to thoseskilled in the art when reading the following detailed description inconnection with the accompanying drawings, wherein:

FIG. 1 is a flow chart illustrating steps of a method in accordance withan embodiment of the invention;

FIG. 2 is a schematic representation of an injector/producerrelationship in a multilayer reservoir having a shale discontinuity,according to an embodiment of the present invention;

FIG. 3 is a schematic illustration of a flow in an injector/producerrelationship illustrating Fermat's principle in determining the locationof shale discontinuities, according to an embodiment of the presentinvention;

FIGS. 4 a-4 c illustrate a testing scenario in a simulated applicationof a method in accordance with an embodiment of the invention;

FIG. 5 a illustrates injection rates, according various embodiments ofthe present invention

FIG. 5 b illustrates measured production rates corresponding to theinjection rates shown in FIG. 5 a;

FIG. 5 c illustrates an estimated FIR and a CM to determine a time delayconstant for the simulated method, according to an embodiment of thepresent invention;

FIG. 6 a illustrates an initial structure, according to an embodiment ofthe present invention;

FIG. 6 b illustrates an output of the algorithm, according to anembodiment of the present invention;

FIG. 6 c illustrates the simulated reservoir structure, according to anembodiment of the present invention;

FIG. 7 shows a comparison between the estimated CM to predict theproduction rate and the measurements, according to an embodiment of thepresent invention;

FIG. 8 a shows the relative well locations for the injector and theproducer wells;

FIG. 8 b shows an estimate probability map for a pilot area, upper side;

FIG. 9 a shows the relative well locations for the injector and producerwells;

FIG. 9 b shows an estimate probability map for a pilot area, lower side;and

FIG. 10 depicts the estimated probability map for the whole field with aknown fault, according to an embodiment of the present invention.

DETAILED DESCRIPTION

In accordance with an embodiment of the present invention, injection andproduction data are used to estimate high permeability channels. Thismay be extended to the 3D case to identify the shale discontinuitieswhich causes cross layer flow. The present method may be implementedwithout stoppage in production schedules.

In particular, the inventors have determined that it may be difficult toidentify shale discontinuities, particular in a layered reservoir withlarge permeability contrast. In order to overcome this difficulty, theinventors have investigated use of flow rates as an indicator. Inparticular, because the shale discontinuities provide channels for thefluid to flow between different layers, when an injection rate ischanged, and producers in the same layer are shut-in, it should bepossible to detect changes in production in different layers. Bycontrolling injection schedule design and including shut-in forproducers, well-pair response time in a 3D layered reservoir may beestimated.

With the injection/production data, a capacitance model (CM) may be usedto characterize the flow property for each well-pair. Specifically, themodel can be viewed as a virtual streamline approach to explain theproduction according to the interaction between well-pairs. In CM, a‘time delay’ constant is introduced to characterize the time delayeffect of the injection signal at the producers. This “time delay”constant is associated with the cross layer travel path, which isdetermined by the location of shale discontinuities. In this regard, theinjector and producer are treated as the transmitter and receiver,respectively, and use the time delay constant as the travel time foreach well-pair. One special property of this problem is that shalediscontinuities may have very different flow properties. This approachmay be referred to as “high contrast travel time tomography.”

After determining a time delay constant, the problem of finding theshale discontinuities can be modeled as a cross-hole travel timetomography. Travel time tomography aims to reconstruct an interiorvelocity model based on measured first-arrival times betweentransmitters and receivers. The velocity model captures the physicalproperties of the region where the signal transmission occurs. Thetransmitted signal can be a seismic wave, an acoustic sound wave, andeven a fluid pressure wave. This technique can characterize large scaleelastic media, in applications such as seismic geophysical explorationand acoustic tomography, in the atmosphere's surface or the oceanlayers. However, the geophysical problem is different from many othertransmission tomographic reconstruction problems (X-Ray CT, PositronEmission Tomography CT, etc.) where the straight-line trajectoryassumption is commonly used. Indeed, in the geophysical problem, theacoustic trajectory may bend severely if the local velocity variation ishigh. Instead of using the straight line trajectory assumption, travelpaths can be characterized by Fermat's principle, where the travel timeobserved corresponds to the path traversed in the shortest time, wherethe time can be obtained as a path integral of the slowness (reciprocalof velocity) function integrated along the travel path.

In many geophysical applications, it is very common to haveheterogeneous structures where the velocity contrast is high, e.g., anarea may have twice the velocity than other areas. Example scenarioswhere this situation can be encountered arise in different applications,e.g., using seismic waves to find permeable fracture zones inground-water flow, monitoring the water/oil saturation in vegetable oilbio-remediation projects, etc. Different from the case that the velocitydistribution only has small variation, the travel path in high velocitycontrast media not only bends severely but is almost dominated by thesehigh velocity structures that form high transmissive channels. Becausethe straight line travel path approximation is not valid, thereconstruction of the velocity model becomes a nonlinear inverseproblem. Furthermore, conventional reconstruction methods are based oniterative linearized algorithms to approximate the travel path. Giventhat in many practical situations, only sparse measured travel time datais available, these methods suffer from problems due to low ray-coverageand severe path bending in the high contrast velocity medium.

Water-flooding is used as a technique in enhanced oil recovery (EOR).Water is injected in a controlled manner in order to provide pressuresupport that can slowly sweep oil into the production wells. During thisprocess, the permeability (measurement of ability to transmit fluid) ofopen fractures can be orders of magnitude higher than that ofsurrounding tight rocks, providing fast pathways for the flow. Thus,travel time through a fracture (which could be modeled as a line in 2Dor a plane in 3D) is much faster than through surrounding rocks. Theflow properties of the reservoir are dominated by these highly‘transmissive’ fracture structures. If a fracture is close to both aninjection well and a production well, most water will directly flowthrough the fracture and thus fail to displace oil in other areas whichis called “water cycling” and thus significantly reduces the efficiencyof oil recovery. Therefore, understanding these fractures is useful inflow characterization to enhance oil recovery.

In one embodiment, high velocity structures (HVS) in a relatively slowhomogeneous background are determined. The velocities are considereddiscrete which calls for solving a ‘discrete’ travel time tomographyproblem. For example, a reservoir can be modeled as a layered structure,where each layer having a different hydraulic conductivity.

Object-based models can be used to represent the HVS. For example,objects can be pre-defined convex polygons, and the geometrical shape ofHVS is represented by a combination of these fundamental objects. Thisallows an arbitrary shape to be represented with convex objects, andprior information about the HVS can be input into the object-basedmodels. For instance, in the fracture characterization problem, thegeometrical shape of fractures can be approximated by lines in 2D modelsor planes in 3D models. The shape of the convex fundamental object canbe defined as a ‘line’. As a result, a small number of objects issufficient to approximate the shape of HVS, which reduces thedimensionality and uncertainty of the problem. Furthermore, the travelpath tracking procedure can be simplified by only considering theshortest time path between objects instead of all grids in the spatialdomain. This approach can reduce the number of unknown variables to theobject parameters and avoids the ‘lack of path coverage’ problem arisingin grid-based models which provides a more stable result in inversion.

In order to estimate the velocity model with the measured travel time, anonlinear inverse problem may be solved. In inverse problems, theestimated velocity model (the solution) may not be unique to the givenmeasurements. One popular approach to handle the non-uniqueness of thesolution is to apply a regularization to favor certain properties in themodel. The regularization methods can lead to solutions that balance thedata-fitting and model-penalty.

In another embodiment, an alternative approach is to estimate theprobability distribution for the model space according to thedata-fitting. This gives a full description of the relativeprobabilities of the different models so that all likely solutions maybe considered.

In one embodiment, the Bayesian approach may be selected to estimate theprobability density map for the model parameter space. Due to the largedimension of the model parameter space, an accelerated random walkalgorithm is provided to explore the model space. In one embodiment,based on the Hamiltonian Monte Carlo method (HMC), an additionalfriction term can be added in the simulation of dynamic system. Thesamples are more likely to fall into local minimum and achieve efficientsampling on high probability regions. The HVS property can then be usedto approximate the neighboring probability distribution. The results arepresented as a probability map showing where the high velocity objectsare more likely to appear in the underlying structure of the system.

In one embodiment, the HVS properties can be used to simplify the pathtracking step, which allows the Hamiltonian Monte Carlo (HMC) samplingprocess to be more efficient. Furthermore, the monotonicity of thetravel time as a function of high velocity object size can be exploitedto approximate the probability distribution in the space of modelparameters.

In one embodiment, the high contrast travel time tomographicreconstruction algorithm described above can be implemented in thegeophysical context and applied in petroleum or gas fields.Understanding the heterogeneous structure of oil and gas reservoirs canprovide information to sweet spot mapping, drilling, productionforecasting and operation efficiency optimization. In enhance oilrecovery (EOR) process, fluids such as water, gas, or chemicals areinjected to increase the amount of oil that can be extracted from thereservoir. The water-flood process is a process where water is used asthe injected fluid. Water injection provides pressure maintenance in thereservoir while displacing part of residual oil to production wells. Tomaximize operation efficiency of the water-flood, steps may be taken touniformly “push” the residual oil where the presence of highpermeability channels affect the horizontal and vertical sweepefficiency. For example, the water-flood process can be applied to atight fractured reservoir. In “tight” fractured reservoir, fracturesprovide conduits for the fluid flow and these conduits provide fastpathways for the injected fluid. Hence, these pathways dominate thesweep efficiency in water-flood.

For example, if a high permeability channel is in close proximity to oneinjection-production well-pair, most injected water may flow throughthis channel and may fail to sweep the oil in other regions. Thisphenomenon is called “water cycling”, which decrease the sweepefficiency and increase the water cut (the ratio of water producedcompared to the volume of total liquids produced). To optimize thewater-flood efficiency, it is desirable to identify the locations ofthese high permeability channels in the reservoir by appropriate changesin injection rates in different wells.

The seismic cross-hole tomography method provides very detailedgeophysical structures. However, in this method, it can be difficult touse the seismic velocity to infer the flow property directly. Othermethods use tracer tests which measure the diffusion process between theinjector-producer pairs and provide the flow permeability distributionin a direct way. A common aspect of these methods is that they allrequire additional equipments and may interrupt daily operations.Seismic cross-hole testing requires inducing seismic waves and deployingsensors to monitor the reflected waves. In tracer tests, chemicals orradioactive materials are injected and concentration response ismonitored to map inter-well flow properties. In addition, repeating thetesting in the same area is complicated and one needs to wait a longtime for the tracer concentration to return back to background levels.

In a water-flood field, the injection/production data is the mostabundant data source and there is freedom to control the injectionrates. Thus, the injection rate can be changed while monitoring thechange in production wells. The estimated response time between wellpairs enables a user to infer the structure of a reservoir using traveltime tomography. In one embodiment, to estimate the travel time, thewater-flood reservoir can be modeled as a linear system by consideringwater injection rates as inputs and total fluid production rates asoutputs. For example, in pulse testing, sensitive differential pressuregauges can be used to monitor the resulting pressure response atadjacent wells.

In another embodiment, the reservoir can be treated as a multi-input andmulti-output system and estimate all the inter-well response by theproduction rates. By using a capacitance model (CM), the relationbetween the injection-production response and the “time delay” constantcan be built. The relation between the injection-production response andthe time delay is approximately proportional to the pressure wavepropagation time. Therefore, the time delay can be estimated from theinjection/production data. The obtained time delay can be used as thetravel time between well pairs. This procedure has the benefit of beingable to be applied without additional cost or significantly affectingdaily operation.

In tight fractured reservoirs, fractures can be viewed as highpermeability channels where the permeability contrast is relativelyhigh, about 10⁵ times permeable than the host rock formation. In thiscase, the travel time for the pressure wave to propagate through thehigh permeability channels can be considered almost negligible. Thetravel time data is restricted to the injection-production well-pairlocations sparsely located in the field. The reconstructed result isrelated to the density of travel ray-path. Thus, the reconstructedreservoir image resolution is fundamentally limited by the spatialdistribution of well locations.

Based on these conditions, the high permeability channels in a reservoircan be characterized by solving a high contrast travel time tomographyproblem with sparse data. If the grid model is used and the conventionaliterative least square methods are applied, the reconstruction resultswould be poor. To identify these high permeability channels anobject-based model which employs “lines” as fundamental objects can beused to represent fractures in 2D.

To estimate the travel time from historical data, a physical model isused to describe the injection/production response. The capacitancemodel can be used to characterize the system response, where the timedelay constant in CM is related to the control volume. When the flowpasses through a homogeneous region, the control volume can be assumedto be proportional to the travel path. By definition, the travel time isthe travel path divided by the local velocity. The time delay constantis proportional to the control volume divided by the productivity index.In a tight fractured reservoir, a simple reservoir model with embeddedhigh permeability channels in a homogeneous background may be used.Thus, if the productivity index in high permeability channels is higherthan the background, the time delay constant inside the highpermeability channels can be ignored. The time delay constant isproportional to the travel time in the high contrast reservoir model.

The capacitance model provides a powerful tool to evaluate thewater-flood performance by treating the injection/production as thesystem input/output and uses signal processing techniques to estimatethe unknown system parameters, which relates to reservoir physicalproperties. The mathematical formulation for one injection-productionwell pair in CM can be expressed by the following equation (1).

$\begin{matrix}{{{\tau \; \frac{q}{t}} + {q(t)}} = {{i(t)} - {\tau*J\; \frac{p_{wf}}{t}}}} & (1)\end{matrix}$

Where i(t) is the injection rate and q(t) is the total production rate.P_(wf) represents the flowing bottom hole pressure (BHP) and J standsfor the productivity index. In the case where the reservoir is matureand the bottom-hole pressure is fixed, then it is possible toapproximate the j-th well production rates based on only thecontribution from injectors by using the following equation (2).

$\begin{matrix}{{q_{j}(t)} = {\sum\limits_{i}{\int_{t_{0}}^{t}{\frac{^{{- {({t - \zeta})}}/\tau_{ij}}}{\tau_{ij}}{I_{i}(\varsigma)}{\varsigma}}}}} & (2)\end{matrix}$

The discrete form of equation (2) can be written as equation (3).

$\begin{matrix}{{q(t)} = {\sum\limits_{k}\left( {\frac{1}{\tau}^{- {({k/\tau})}}{i\left\lbrack {t - k} \right\rbrack}} \right)}} & (3)\end{matrix}$

Where t is equal to 1, 2, . . . , N. For multiple injection-productionwell pairs, additional parameters l_(ij) for inter-well connectivity,which represent flow distribution from injector i to different producersj. The j-th well production rate is the sum of flow from differentinjector wells i. This can be expressed by the following equation (4).

$\begin{matrix}{{q_{j}(t)} = {\sum\limits_{i}{\sum\limits_{k}{l_{ij}\left( {\frac{1}{\tau}^{- {({k/\tau})}}{i\left\lbrack {t - k} \right\rbrack}} \right)}}}} & (4)\end{matrix}$

By using equation (4), a linear multi-input multi-output (MIMO) systemcan be used to model the relation between injection/production rate. Theimpulse response is controlled by the “time delay constant” τ_(ij)between each injector-producer pair which is defined by the totalcompressibility c_(t), the productivity index J, and the pore volume Vpwhich are associated with the control volume between injector i andproducer j pair. The time delay constant τ_(ij) between eachinjector—producer pair can be expressed by the following equation (5).

$\begin{matrix}{\tau_{ij} = \left( \frac{c_{t}V_{p}}{J} \right)_{ij}} & (5)\end{matrix}$

From above, it is clear that if there is an area with higherproductivity index J, the time delay constant will be much smaller. Thisimplies that if the time delay constant for a specific well-pair is verysmall, it is probable that a majority of the flow path is through a highpermeability layer.

For a homogeneous reservoir, the control volume is roughly proportionalto the distance between the injector and producer. When a cross-layerflow in a layered reservoir is considered, the travel path will be thepath with the least time cost in accordance with Fermat's principle.This implies that the flow will take advantage of the higherpermeability layer.

A total control volume can be modeled as the cascade of several smallcontrol volumes, e.g., for example two control volumes which can be twodifferent layers. Each control volume is proportional to a line segmentwhich represents the flow path through the control volume. Therefore,the time-delay constant between a well pair will be the summation of twoparts: the time delay constant of the first layer (e.g., upper layer)and the time delay constant of the second layer (e.g., lower layer).Each is given by the ratio of flow path in different layers, which isdetermined by the location of the shale discontinuities because theyprovide the entries for the fluid to travel through different layers.

In order to get a reliable estimate of the parameters for CM, aPseudo-Noise (PN) sequence may be used as the injection schedule. Byusing PN sequence as the input, it may be possible to achieve the lowestcovariance for the parameter estimation. In addition, by using PNsequences as inputs, lowest average cross-correlation between inputs mayalso be achieved.

When the injection rate in the lower permeability layer is changed,usually it is possible to detect the change of production even for theproducers located in the higher permeability layer. However, in someinstances, it may be difficult to detect this cross-layer flow if theinjection rate in the high permeability layers is changed. This isbecause the high permeability layer provides a fast pathway for thefluid so only a small amount will flow through the shale discontinuityto the low permeability layer. Therefore, when designing the injectionschedule it may be useful to shut-in the producers in the same layer toforce the fluid to flow through different layers. In this situation, thecross layer flow can be detected by monitoring the change of production.

The algorithm to identify the shale discontinuities with cross-layertime delay constant can be referred as the “high contrast travel timetomographic reconstruction”. One of the issues with this approach isthat the sensor locations may be limited in the injectors/producers, andthe flow velocity contrast of the shale discontinuities may be veryhigh. If the conventional grid based inversion method is used, theresult may result in a blurred image and it may be difficult to identifythe shale discontinuities.

As such, instead of using the grid based method, an object-basedrepresentation of the shale discontinuities may be chosen. Thereconstruction algorithm will modify the object parameters for the shalediscontinuities to match the measured delay time constant. It can bedescribed as a ‘Forward-Backward’ iterative procedure.

FIG. 1 is a flow chart illustrating steps of a method for determiningthe time delay constant, according to an embodiment of the presentinvention. In the forward step, a time delay is calculated based on thecurrent or initial locations of shale discontinuities based on a currentor initial reservoir model, at S10. At S12, a test is performed todetermine whether or not a difference between the time delay and ameasured time delay is small enough, i.e., smaller than a set threshold(e.g. smaller than 1% to 5% of the measured value). If the calculatedresult fits the measurement (i.e., the difference between themeasurement and the current estimate is small or smaller than the setthreshold), then the selected current reservoir model can be consideredas a solution and the algorithm is stopped, at S14. If not, the‘backward’ step is performed to update the locations of shalediscontinuities (i.e., the initial reservoir model is updated) based onthe difference of the current predicted response time with themeasurements, at S16. Therefore, if the difference is greater than theset threshold, the method includes updating the initial reservoir modelbased on the difference, calculating the time delay based on the updatedreservoir model; and determining whether or not a difference between thetime delay based on the updated reservoir model and the measured timedelay is smaller than the threshold. This loop procedure can be repeateduntil convergence between the calculated time delay and the measuredtime delay is achieved. That is, the loop procedure is repeated untilthe difference between the calculated time delay and the measured timedelay is smaller than the set threshold.

The detail of the “forward” step can be described as follows. With thecurrent location of shale discontinuities, the flow travel path can becalculated by Fermat's principle. The shale discontinuities act as entrypoints to the different layers, and the time delay constant can beapproximated by the cascade of different layers which is proportional tothe travel distance inside. In this case, the time delay constant iswell defined by the locations of shale discontinuities.

FIG. 2 is a schematic representation of an injector/producerrelationship in a multilayer reservoir having a shale discontinuity,according to an embodiment of the present invention. Injector well 10and producer well 12 provided within rock formation 13 are separated bydistance L. As illustrated in FIG. 2, there two shale layers 14A and14B. Shale layer 14A has a lower permeability and shale layer 14B has ahigher permeability. As shown in FIG. 2, there is a shale discontinuity16 between shale layer 14A and shale layer 14B, the discontinuitycorresponding to a fracture in the rock formation 13 where fluid flow(indicated schematically by the arrows) is privileged. For example, ifit is assumed that the flow velocity of the layer (shale layer 14A)having lower permeability is equal to 1 (VL=1) and the flow velocity ofthe layer (shale layer 14B) having a higher permeability is 10 (VH=10),the flow path can be separated as two segments belonging to the twodifferent layers 14A and 14B and the delay time can be calculated usingthe following equation (6).

Travel time=L ₁/1+L ₂/10, with L ₁ +L ₂ =L  (6)

where L represents the geometrical distance between theinjector-producer well 10 and 12, and the ratio of two segments iscontrolled by the location of shale discontinuities. The distance L1corresponds to the distance traversed by the injected fluid throughlayer 14A. The distance L2 corresponds to the distance traversed by theinjected fluid through layer 14B. The travel time within thediscontinuity or fracture 16 is considered small relative to the traveltimes of fluid flow within the layers 14A and 14B. Therefore, the traveltime within fracture 16 is not taken into account when calculating thetotal travel time using equation (6).

FIG. 3 is a schematic illustration of a flow in an injector well andproducer well relationship illustrating Fermat's principle indetermining the location of shale discontinuities, according to anembodiment of the present invention. Injector well 20 and producer well22 provided within rock formation 23 are separated by distance L. Asillustrated in FIG. 3, there are two types of shales 24A and 24B. Shale24A has a lower permeability and shale 24B has a higher permeability. Asshown in FIG. 3, there is a shale discontinuity 26 between shale 24A andshale 24B, the discontinuity corresponding to a fracture in the rockformation 23. In this model of the earth, the relationship between thevarious shales (i.e., shales with various permeabilities) and theposition of the discontinuity is not as straight forward as in theexample shown in FIG. 2. Therefore, in order to determine an appropriatemodel that provides an estimate of the travel time and thus the locationof the discontinuity, the backward step in the loop S16 is performed. Inthe ‘backward’ step, the locations for the shale discontinuities 26 isupdated base on the mismatch of the measured and predicted time delayconstant τ. Time delay constant τ is the sum of time delay τ(1) withinshale 24A and time delay τ(2) within shale 24B. Time delay τ(1) isproportional to the traversed length L1 within shale 24A and inverselyproportional to the flow velocity within shale 24A. Similarly, timedelay τ(2) is proportional to the traversed length L2 within shale 24Band inversely proportional to the flow velocity within shale 24B.

Because the shale discontinuities are represented by objects withparameters, we can formulate the updating as a parameter estimationproblem. Current object parameter can be expressed by the followingequation (7).

$\begin{matrix}{\theta = {\underset{\theta}{Argmin}{{{T(\theta)} - t}}}} & (7)\end{matrix}$

where T(θ) is the calculated response time based on current objectparameters θ and t is the measurement. A gradient search is used to findthe optimal parameters to represent the location of the shalediscontinuities.

FIGS. 4 a-4 c illustrate a testing scenario in a simulated applicationof a method in accordance with an embodiment of the invention. A twolayer reservoir model with high permeability contrast was built and acommercial simulator (CMG) is used to test the method. The testing casewas a line drive with 5 injectors and 5 producers. Three of theinjectors/producers are in the low permeability layers, and the twoother injectors/producers are in the high permeability layers. A simpleshale discontinuity was placed and high permeability assigned to it tosimulate the effect. For the reservoir model configuration, a highpermeability channel is simulated. In the simulation, the highpermeability channel connects injector well 5 with producer well 3. Apulse testing is applied in the simulation. The pulse testing uses astep or square function or wave as the change of injection rate. Theinjection rate is changed each time so the response can be easilyestimated from the change of production rate. The travel time is definedby the time to reach 90% of the final change of the production rate.FIG. 4 a shows the configuration of injector wells I2 and I4 relative toproduction wells P2 and P4 in the first layer. FIG. 4 b shows theconfiguration of injector wells I1, I3 and I5 relative to producer wellsP1, P3 and P5 in the second layer. FIG. 4 c illustrates the shale layersand their respective permeabilities (shown on the scale strip) withinthe rock formation and the relative position of the 5 injector wells I1through I5 and 5 producer wells P1 through P5.

A PN sequence was used as the injection schedule and the producers inthe same layer were shut-in for the testing period. The changes ofproduction can be measured for the producers located in differentlayers. A time delay constant was retrieved in CM byinjection/production data, and used to identify the location of shalediscontinuities as shown in FIGS. 5 a-5 c.

FIG. 5 a depicts two plots with two different injection rates. FIG. 5 bdepicts two plots with two production rates corresponding to the twoinjection rates shown in FIG. 5 a. FIG. 5 c represents a plot of anestimated FIR and a CM to determine a time delay constant for thesimulated method. The plot represents a correlation between injectionand production as a function of time (the horizontal or x-axisrepresenting time and the vertical or y-axis representing the normalizedresponse). As it can be noted, the average time delay to obtain a changeof production at the production site is in this case about 5 to 6 days.

To apply the method in a real field, one of the useful parameters is thesampling frequency. Another useful parameter is accuracy of the data.Usually it is possible to obtain reliable and frequent injection data,but production rates are obtained from the well-test, which canperformed periodically, e.g., on a weekly basis. This factor tends tolimit the time-resolution when the time delay constant is measured. Forexample, assuming the time delay constants for the well-pairs are 1 and4 days, when injection rates are increased, it is possible to measurealmost the same changes in production rate if it is measured bywell-test. To deal with this problem, it is possible to use an inferredproduction value provided by a controller pump which provides morefrequent production data.

FIG. 6 a illustrates an assumed initial structure of the highpermeability channel or discontinuity, according to an embodiment of thepresent invention. The assumed initial structure is used in the initialmodel to estimate an initial time delay. FIG. 6 b illustrates an outputof the algorithm, according to an embodiment of the present invention.FIG. 6 b depicts the position of the high permeability channel or shalediscontinuity obtained through refining the initial model until thedifference between the estimate time delay and the measured time delayis relatively small and below a certain threshold. FIG. 6 c illustratesthe simulated reservoir structure, according to an embodiment of thepresent invention. The scale strip next to the graph provides a scale ofthe permeability of the rock formation. The bar within the graphindicates the position of the discontinuity. A good match is obtainedbetween the measurement and the calculated position or location of thediscontinuity (fault) as well as the direction of the discontinuity(fault).

A field experiment under water-flood in an oil field is performed. Thefield is considered to be a tight fractured reservoir. In theexperiment, a pilot area with 12 injection wells is considered. The wellinstallation is a line-drive and roughly parallel to the direction offractures which are estimated from seismic survey.

Two types of systems in the oil field are used to provide the wellproduction data. One is the “well-test” data, which accumulates a threephase production rate of one specific well for a short period of time.However, many wells share the same separation facility. Therefore, thesewells are tested sequentially. This limits the sampling rate of theproduction rates, because the sampling period is the sum of the testingperiod for all production wells. In this case, one measurement pointevery two weeks is taken which makes our production rate sparselysampled in the well-test data.

The other data source is the pump-of-control (POC) data, which measuresthe pump load change in every stroke. If the number of strokes countedis multiplied by the load change, the theoretical total daily productioncan be calculated. To get a reliable data with high sampling rate, thePOC data is used and the pumping efficiency is calibrated with thewell-test data. After calibrating with the pumping efficiency, the POCdata would be able to provide a daily production rate. However, we doobserve unusual low production rate, which might be outliers due to thepower outage or loss of transmission. Thus, in the linear modelestimation we choose l₁ instead using l₂ norm as the error metric, whichis more balance between fitting good data points and outliers.

To estimate the CM parameters, the direct approach is chosen in thiscase. The reason is because the input data length is limited (about 90),and thus, the multi-stage approach is not suitable because the unknownparameters in FIR model will be much larger than the measured data. Froma previous discussion in the above paragraphs, it is established thatthe direct estimation uses fewer unknown variables but needs to solve anonlinear optimization. Therefore, the estimated CM might be a solutionin a local minimum and not represent the system well. To verify theestimated CM, cross-validation is used which divides the measured datainto two parts. The first part is the training data, which is used toestimate the CM parameters. The second part is the testing data, whereestimated CM is used to predict the production rate and compare with themeasurements. If the error is too high, the estimated CM is determinedto be not accurate and the non-linear optimization is re-run. FIG. 7shows a comparison between the estimated CM to predict the productionrate and the measurements, according to an embodiment of the presentinvention.

Following the same principle in numerical simulation, the travel time isonly used when the inter-well connectivity is greater than 10%. Thepilot area is separated into upper and lower side and the travel timebetween well pairs inside each region is estimated. The reason why theestimation for the well pair is not performed for the whole field isthat that solving for the whole field needs a greater number ofvariables. In the line drive water-flood, it can be assumed the flowonly comes from the nearby two rows of injectors. Hence, solving twosmaller fields is selected instead.

FIG. 8 a shows the relative well locations for the injector and theproducer wells. FIG. 8 b shows an estimate probability map for pilotarea 1, upper side. It can be noted that the estimated high permeabilityis roughly parallel to the installation of the wells which is consistentwith the prior information of fracture direction. FIG. 9 a shows therelative well locations for the injector and producer wells. FIG. 9 bshows an estimate probability map for pilot area 1, lower side. Similarto FIG. 8 b, it can also be noted that the estimated high permeabilityis roughly parallel to the installation of the wells which is consistentwith the prior information of fracture direction.

The estimate travel time results for the upper part of pilot area isshown in FIG. 8 b, and the lower part is shown in FIG. 9 b. In the upperpart, two reliable production rate measurements are obtained, and thetravel time used are t_(1,1), t_(3, 2), t_(4, 2), t_(5,1), t (5, 2),t_(6,1), t_(7,1), t_(7,2), t_(8,1), t_(8,2), a total of 8 well pairs.The result shows in the upper side, the estimated high permeabilitychannels are close to production well 2 and roughly parallel to theinstallation of wells. This is consistent with the fact that the wellproduction is a 24 hours run, high production well, which implies thatthere are high permeability channels nearby. In the lower part, oneproduction well is present and the estimated travel time are t_(7,3),t_(9,3), t_(10,3), t_(11,3), and t_(12,3) where there are fewer measuredtravel times to estimate the high permeability map. The result in twoareas are then combined and compared with the known fault location. Theresult shows the appearance probability of high permeability channels inthe fault zones is very low, and roughly parallel to the wellinstallation which agrees with the seismic survey. FIG. 10 depicts theestimated probability map for the whole field with a known fault,according to an embodiment of the present invention. The position of theknown fault is indicated on the probability map. The results show thepermeability channel has a very low appearance probability in the faultzone consistent with the prior survey.

The “time delay constant” can be used in a CM method to detect highpermeability channels. The method uses the injection-production. Themethod uses the change of injection as the active probing and detect thefield changes in real-time without altering the average dailyproduction. In order to apply this to real field data, some practicalissues are noted. First, the data sampling period and quality can playan important role. Usually, a reliable daily injection rate data isobtained. However, production rates are often obtained from bi-weeklywell-test data. This may reduce the time resolution of estimated traveltimes. The POC data calibrated with well-test data can be used to getthe daily production rate. In addition, in practical applications, thedistribution of well locations is related to spatial resolution. Forexample, in an extreme case where one injector-producer pair is in thehorizontal direction, any high permeability channel exactly in thevertical direction will not affect the lag time. Therefore, the highpermeability channels may not be visible in this situation. In order todetect the high permeability channel in any arbitrary direction, thewells can be uniformly located in the field and to cover all angles.

As will be appreciated, the method as described herein may be performedusing a computing system having machine executable instructions storedon a tangible medium. The instructions are executable to perform eachportion of the method, either autonomously, or with the assistance ofinput from an operator. In an embodiment, the system includes structuresfor allowing input and output of data, and a display that is configuredand arranged to display the intermediate and/or final products of theprocess steps. A method in accordance with an embodiment may include anautomated selection of a location for exploitation and/or exploratorydrilling for hydrocarbon resources. Where the term processor is used, itshould be understood to be applicable to multi-processor systems and/ordistributed computing systems.

Those skilled in the art will appreciate that the disclosed embodimentsdescribed herein are by way of example only, and that numerousvariations will exist. The invention is limited only by the claims,which encompass the embodiments described herein as well as variantsapparent to those skilled in the art. In addition, it should beappreciated that structural features or method steps shown or describedin any one embodiment herein can be used in other embodiments as well.

1. A method of locating discontinuities in a reservoir, comprising:identifying a plurality of injection and production sites in thereservoir; injecting fluid into selected injection sites in accordancewith an injection schedule; monitoring an output at the productionsites; determining a time delay between the selected injection sites andthe production sites based on the monitored output; and determining alocation of discontinuities based on the determined time delay.
 2. Themethod according to claim 1, wherein determining the time delaycomprises determining the time delay using a forward/backward iterativeprocedure.
 3. The method according to claim 1, wherein determining thelocation of discontinuities comprises using a tomographic method and thedetermined time delay.
 4. The method according to claim 3, wherein usingthe tomographic method includes reconstructing an interior velocitymodel with a rock formation based on measured arrival times between atransmitter and a receiver.
 5. The method according to claim 4, whereinthe velocity model captures physical properties of a region of the rockformation where a transmitted signal propagates, the transmitted signalincluding a seismic wave, an acoustic wave, or a fluid pressure wave. 6.The method according to claim 1, further comprising closing selectedproduction sites that are in a common layer with the selected injectionsites to force the injected fluid to flow between different layers. 7.The method according to claim 1, wherein determining the time delaybetween the injection sites and the production sites based on themonitored output comprises: calculating a first time delay based oninitial locations of discontinuities that are based on initial reservoirmodel; determining whether or not a difference between the calculatedfirst time delay and a measured time delay is smaller than a setthreshold; if the difference is smaller than the set threshold, then theinitial reservoir model is used to determine the locations of thediscontinuities.
 8. The method according to claim 7, wherein determiningthe time delay further comprises: if the difference is greater than theset threshold, updating the initial reservoir model based on thedifference; calculating a second time delay based on the updatedreservoir model; and determining whether or not a difference between thesecond time delay based on the updated reservoir model and the measuredtime delay is smaller than the threshold.
 9. The method according toclaim 8, wherein the updating of reservoir model, the calculating of thesecond time delay and the determining whether or not the differencebetween the calculated second time delay based on the updated reservoirmodel is repeated until the difference between the calculated secondtime delay and the measured time delay is smaller than the setthreshold.
 10. The method according to claim 7, the set threshold isequal to approximately one percent to 5 percent of the measured timedelay.
 11. A system for locating discontinuities in a reservoir,comprising: a processor configured to: determine a time delay betweeninjection sites and production sites in a the reservoir based onmonitored output at the production sites from injecting fluid toselected injection sites in accordance with an injection schedule; anddetermine a location of discontinuities based on the determined timedelay.
 12. The system according to claim 11, wherein the processor isconfigured to determine the time delay using forward/backward iterativeprocedure.
 13. The system according to claim 11, wherein the processoris configured to determine the location of discontinuities using atomographic method and the determined time delay.
 14. The systemaccording to claim 13, wherein the tomographic method includesreconstructing an interior velocity model with a rock formation based onmeasured arrival times between a transmitter and a receiver.
 15. Thesystem according to claim 14, wherein the velocity model capturesphysical properties of a region of the rock formation where atransmitted signal propagates, the transmitted signal including aseismic wave, an acoustic wave, or a fluid pressure wave.
 16. The systemaccording to claim 1, wherein the processor is configured to: calculatea first time delay based on initial locations of discontinuities thatare based on initial reservoir model; determine whether or not adifference between the calculated first time delay and a measured timedelay is smaller than a set threshold; if the difference is smaller thanthe set threshold, then the initial reservoir model is used to determinethe locations of the discontinuities.
 17. The system according to claim16, wherein the processor is configured to: update the initial reservoirmodel based on the difference, if the difference is greater than the setthreshold; calculate a second time delay based on the updated reservoirmodel; and determine whether or not a difference between the second timedelay based on the updated reservoir model and the measured time delayis smaller than the threshold.
 18. The system according to claim 17,wherein the processor is configured to repeat the updating of reservoirmodel, the calculating of the second time delay and the determiningwhether or not the difference between the calculated second time delaybased on the updated reservoir model until the difference between thecalculated second time delay and the measured time delay is smaller thanthe set threshold.
 19. The system according to claim 16, the setthreshold is equal to approximately one percent to 5 percent of themeasured time delay.